Hydrogen diffusion in the proton conductor Gd-doped barium cerate 



Jessica Hermet 1 ' 2 , Marc Torrent 1 , Frangois Bottin 1 , Guilhem Dezanneau 2 , and Gregory Geneste^ 

1 CEA, DAM, DIF, F-91297 Arpajon, France 
2 Laboratoire Structures, 
Proprietes et Modelisation des Solides, 
UMR CNRS 8580, Ecole Centrale Pans, 

Grande Voie des Vignes, 
92295 Chdtenay-Malabry Cedex, France 
(Dated: January 3, 2013) 

The energy landscape and diffusion barriers ol protonic defects in Gd-doped BaCe03, a compound 
candidate as electrolyte for protonic ceramic fuel cells, have been investigated by density functional 
theory calculations, starting from a previously computed energy landscape consisting of 16 kinds 
of stable sites (8 close to dopants and 8 far from them). The simplified string method has been 
used to determine accurately the Minimum Energy Paths between those sites, that might imply 
either proton reorientations, intra-octahedral or inter-octahedral hopping mechanisms. At contrast 
with simple cubic perovskites such as barium stannate or barium zirconate, very different values for 
energy barriers (from 0.02 eV to 0.58 eV) are found in this highly distorted orthorhombic perovskite, 
and no specific process appears to be clearly rate-limiting. Some inter-octahedral hoppings (when 
possible) are found to be more favourable than the intra-octahedral ones, while reorientations exhibit 
a wide range of energy barriers. 



PACS numbers: 

I. INTRODUCTION 

Since the discovery of protonic conductivity in 
aliovalent-doped SrCeO^^I, protonic conduction in 
perovskite-type oxides ABO3 has been the subject of 
numerous studies, experimental as well as computa- 
tional^. The high protonic conductivity in perovskite 
oxides opens the way for a wide range of technological ap- 
plications such as Protonic Ceramic Fuel Cells (PCFCs), 
hydrogen separators, etc. However, if the diffusion of 
protons has been extensively explored by ab initio calcu- 
lations in cubic perovskites such as barium zirconate^', 
only very few works have studied this phenomenon in or- 
thorhombic perovskites^, although excellent proton con- 
ductors, such as SrCe03 or BaCe03, can be found among 
such systems. 

Proton conductors are usually obtained by replacing 
some cations of a host oxide compound by cations with 
lower valence. In perovskite oxides having a tetravalent 
element on the B site (Ti, Zr, Ce, Sn), this can be done 
by inserting on this site a trivalent element. Such sub- 
stitution creates charge-compensating oxygen vacancies 
that make the compound reactive with respect to wa- 
ter dissociation if it is put in contact with humid atmo- 
sphere. Such hydration reaction is commonly written, 
using Kroger- Vink notations, as 

H 2 + VS' + Oo -> 20H o . (1) 

It generates protonic defects OHq, localized approxi- 
mately along [100]-type directions inside the interoctahe- 
dral space of the perovskite network, and that can move 
from an oxygen site to another by simple thermal ac- 
tivation. Three possible motions of the proton in the 
perovskite network have been distinguished: 



(i) the reorientation: the OH bond does not break and 
simply turns by 90° around the B-O-B axis containing 
the oxygen atom. 

(ii) the intra-octahedral hopping: the proton leaves its 
oxygen site to move on another oxygen site of the same 
octahedron. 

(iii) the inter-octahedral hopping: the proton leaves its 
oxygen site to move on another oxygen site that does not 
belong to the same octahedron. 

In a previous worM^l, we have studied by density- 
functional theory calculations the thermodynamics of 
hydration and oxidation of Gd-doped barium cerate 
BaCei_,5Gd,50 3 _£ (BCGO). In particular, we have 
showed that hydration was an exothermic process and 
accurately determined the energy landscape of the pro- 
ton near and far from the Gd dopant. We showed that 
this energy landscape can be well approximated by a sur- 
face with 16 kinds of local minima (8 in the close vicin- 
ity of the dopant, and 8 further). This complexity is 
the consequence of the highly distorted geometry of the 
host BaCeOs, that adopts in its ground state the Pnma 
space group. Consequently, proton migration through- 
out such energy surface involves many different energy 
barriers that need to be explored in order to get insight 
into proton conduction at the macroscopic scale. Pre- 
vious works have studied proton migration in BaCeOs, 
but only in the cubic phase^OU. Therefore, in this work, 
we present an exhaustive study of the Minimum Energy 
Paths associated to the possible motions for the proton 
in orthorhombic BCGO, and the values of their energy 
barriers. 
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II. COMPUTATIONAL DETAILS 
A. Method 

We have performed density functional the ory (D FT) 
calculations using the plane- wave code ABINIT 15 16 . The 
Generalized Gradient Approximation (GGA-PBE-^l) was 
employed to describe electronic exchange and correlation. 
The calculations were carried out in the fram ework of the 
projector augmented wave (PAW) approacfPEH. The 
same supercell as that of Ref. [TT] was used: it consists 
of 80 atoms and has an orthorhombic symmetry (Puma 
space group). The First Brillouin Zone of this supercell 
was sampled by a 2x2x2 k-point grid, and the plane- 
wave cutoff was set to 20 Ha. The numerical accuracy 
on the total energies associated to this scheme is better 
than 1 mHa/atom. The cut-off radii of our PAW atomic 
data can be found in Ref. [TT1 

In order to compute Minimum Energy Paths, the first 
task was to identify the stable sites of the proton in 
BCGO, which was previously achieved in Ref. [Tl . This 
was performed by substituting in the 80-atom super- 
cell one Ce by one Gd and introducing one hydrogen 
atom, that was placed in its different possible sites, 
close to the Gd dopant and far from it. In each con- 
figuration, the atomic positions were optimized until all 
the cartesian components of atomic forces were below 
lxlO" 4 Ha/Bohr (« 0.005 eV/A). 

The possible energy barriers between pairs of stable 
protonic sites have then bee n co mputed using the so- 
called simplified string method^lSU. The simplified string 
method is an iterative algorithm allowing to find the Min- 
imum Energy Path (MEP) between two stable configu- 
rations. It consists in discretizing the path into equidis- 
tant configurations, that we call "images". At each iter- 
ation, a two-step procedure is applied: first, each image 
is moved along the direction given by the atomic forces 
(evolution step) , then the images are redistributed along 
the path in order to be kept equidistant (reparametriza- 
tion step). To determine the number of iterations of 
string method, we used an optimization criterion related 
to the energy of the images: the optimization of the 
MEP is stopped when the total energy (averaged over 
all the images) difference between an iteration and the 
previous one is lower than lxlO -5 Ha. In such an al- 
gorithm, the result should be carefully converged with 
respect to the number of images along the path, which 
forced us to use up to 19 images in the case of some 
intra-octahedral hopping processes. Once the MEP has 
been correctly converged, the maximum energy along the 
path provides us the transition state, and thus the energy 
barrier of the corresponding process (hopping or reorien- 
tation). Finally, we point out that all the atoms of the 
supercell were allowed to move during the computation 
of the MEP, thus providing energy barriers in a "fully- 
relaxed" system. 

For the sake of numerical efficiency, we have used the 
three traditional levels of parallelization present in the 



ABINIT code (k-points, bands, plane waves) together 
with a fourth level on the images of the system used to 
discretize the MEP. This fourth level has a quasi-linear 
scalability and, since the number of images used to dis- 
cretize the path can be as large as 19, thousands of cpu 
cores can be used to compute and relax the MEP with 
high efficiency. Typical jobs were done on 3000 cpu cores 
using these four parallelization levels, allowing us to take 
maximal benefit of the potentialities of parallel super- 
computers. 



B. Approximations to the computation of energy 
barriers 

Additional remarks have to be mentionned about the 
limitations of our approach and the approximations made 
to compute the energy barriers. 

First of all, the string method, like the Nudged Elastic 
Band method, allows to compute the Minimum Energy 
Paths between two stable configurations and thus to ob- 
tain "fully-relaxed" (static) barriers, as opposed to "dy- 
namical" barriers that would be obtained, for instance, 
by counting the occurences of each event within a molec- 
ular dynamics run and fitting the rates by an Arrhenius 
law. Static barriers neglect some collective effects and 
the so-called recrossing processes. In theory, they make 
sense only if the whole structure is able to relax instan- 
taneously when the proton moves from a stable position 
to another. However, the time scale associated to the 
hydrogen motion is much smaller than the ones of the 
deformation of the surrounding structure, which involves 
much softer phonon modes. The motion of protons in an 
unrelaxed envionment would naturally lead to higher bar- 
riers than those calculated from fully-relaxed DFT calcu- 
lations. Nevertheless, as shown by Li and Wahnstrom 22 
in metallic palladium, the jump of the proton has to be 
considered in a reverse way. Due to the vibrations of 
surrounding atoms and to the high vibration frequency 
of hydrogen, protons currently jump at a moment where 
the surrounding atoms are in a geometrical configuration 
close to the calculated relaxed one. That is why the cal- 
culated barriers can result very close to those currently 
observed. Further work should be nervertheless neces- 
sary to verify that the proton jump, for instance during 
ab inito Molecular Dynamics simulations, occurs for a 
geometry of surrounding atoms close to that calculated 
in the fully relaxed DFT static scheme. 

Second, the present barriers do not include quantum 
contributions from zero-point motions. They are valid in 
the limit where nuclei can be considered as classical par- 
ticles. If this approximation is correct for heavy atoms in 
the temperature range interesting PCFCs, this is not so 
obvious for the protorP. Indeed, proton tunnelling might 
occur and thus significantly lower the barrier height, es- 
pecially in the hopping case^. This approximation leads 
to overestimated barriers. 

Last, the use of the Generalized Gradient Approxima- 
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tion tends to underestimate the activation energy for pro- 
ton transfer in hydrogen-bonded systems^. This under- 
estimation is due to an over-stabilization of structures in 
which an hydrogen is equally shared between two elec- 
tronegative atomiP^. 

Consequently, the barriers presented in this work 
purely reflect the GGA potential energy surface of the 
proton in its host compound. They are static barriers, 
free from collective, dynamical and quantum effects. 



III. REVIEW OF PRELIMINARY RESULTS: 
STRUCTURE OF BACEO3 AND PROTONIC 
SITES 

A. BaCe0 3 and BCGO structure 

As many perovskites^, BaCe03 has an orthorhombic 
structure (Pnma space group^) at room temperature 
(RT). At high temperature, it undergoes three structural 
phase transitions, the first one at ss 550 K towards an 
Imma structure, and the second one at ss 670 K towards 
a rhombohedral R3c structure. At very high tempera- 
ture (~ 1170 K), it eventually evolves towards the parent 
PmZm cubic structure, that of the ideal perovskitc. The 
presence of dopants randomly distributed throughout the 
matrix may change transition temperatures. However 
Melekh et alW^ found that the first transition in 10%-Gd- 
doped BaCeOs occurs around 480-540 K, close to the one 
they found for pure BaCe0 3 of 533 K. At RT, Gd-doped 
BaCeOs is therefore orthorhombic. 

Our calculations provide optimized configurations and 
Minimum Energy Paths. These computations are thus 
relevant when performed in combination with the ground 
state structure of BaCeOs, i.e. the orthorhombic Puma 
structure, which was used as starting point in all the cal- 
culations, and globally preserved along the optimizations 
procedures. The computed energy barriers can therefore 
be used to understand proton diffusion in BCGO below 
~ 550 K. However, from a more general point of view, 
the present results provide a useful microscopic insight 
into proton diffusion in a low-symmetry perovskite com- 
pound, typical of those used as electrolytes in Proton 
Ceramic Fuel Cells (the Puma structure is common to 
many perovskites such as cerates, zirconates, titanates or 
stannates) . 

The structural parameters obtained for BaCe03 and 
BCGO within the present scheme can be found in Ref.fTTl 
They are in excellent agreement with experiments, de- 
spite a slight overestimation of the lattice constants re- 
lated to the use of the GGA. 



B. Protonic sites in perovskites: general 
considerations 

As previously explained, proton conduction in an 
ABO3 perovskite compound - where B is a tetravalent 



element - might be obtained by substituting B atoms 
by trivalent elements such as Gd (this creates oxygen 
vacancies by charge compensation) and by subsequently 
exposing the new compound to humid atmosphere. The 
protons as charge carriers then appear through the disso- 
ciation of water molecules into the oxygen vacancies, ac- 
cording to the well-known hydration reaction (see equa- 
tion [l]). 

The precise location of the stable protonic sites in the 
perovskite network seems to strongly depend on the lat- 
tice parameter and distortion of the host compound. It 
is commonly admitted that protons are bonded to an 
oxygen atom and remain in the form of hydroxyl groups 
located on oxygen sites. But the orientation of the O-H 
bond is not that clear. On the one hand, it was proposed 
that it could be oriented alon g the BOg octahedra edge 
because of its dipolar momentP^l leading to 8 possible 
sites per oxygen atom. On the other hand, previous ex- 
perimental and ab initicE^U studies have found only 
four sites per oxygen atom oriented along the pseudo- 
cubic directions. 

In fact, the stable protonic sites seem to be indeed 

(i) along or close to the octahedra edge for per- 
ovskites with relatively small lattice constant do, such 
as SrTiO;^ 1 or LaMnO^ (o =3.91 A), leading to the 
existence of 8 protonic sites per oxygen atom, 

(ii) along the pseudo-cubic directions for perovskites 
with l arge lat tice constant, such as SrZrOj^ or 
BaCeOj 1 - 1 1 30 1 32 1 (pseudo-cubic lattice constant a =4.14 
and 4.41 A respectively), leading to the existence of 4 
protonic sites per oxygen atom. 

This trend can easily be explained : as the lattice con- 
stant decreases, the nearest oxygen gets closer and closer 
to the proton, attracting it sufficiently (through hydro- 
gen bond) to bend the O-H bond towards the octahedron 
edge. 



C. Protonic sites in Gd-doped BaCeOa 

In our previous calculations on BCGCW, which has a 
large pseudo-cubic lattice constant of 4.41 A, we found 
indeed four stable protonic sites per oxygen atom. Con- 
sidering that the Pnma structure contains two inequiva- 
lent oxygen atoms 01 and 02, this leads to the existence 
of 8 inequivalent stable positions for the proton, if we 
ignore the symmetry-breaking caused by the presence of 
dopants. These positions have been labeled la, lb, lc 
and Id for those attached to 01 (apical oxygen), and 
2a, 2b, 2c and 2d for those attached to 02 (equatorial 
oxygen), see Fig. [I] 

However, when one Ce atom is replaced by a Gd 
dopant, both the translational symmetry and the sym- 
metry between the four equatorial oxygens 02 of the 
first coordination shell of this specific B site are broken. 
More precisely, the presence of Gd splits the four 02 
into two pairs of symmetry equivalent oxygen atoms (02 
and 02'). The four inequivalent protonic sites related 
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FIG. 1: The 8 stable positions for the proton around the Gd 
dopant. 



to 02 (2a, 2b, 2c and 2d) are thus split into 8 inequiv- 
alent sites, called 2a, 2b, 2c, 2d, 2a', 2b', 2c' and 2d'. 
The first coordination shell of Gd exhibits therefore 12 
kinds of inequivalent protonic sites. Beyond this shell, 
the symmetry-breaking is even more complex. 

Nevertheless, we have shown in Ref. [TT| that this new 
emerging complex protonic energy landscape can be very 
well approximated by a surface containing 16 kinds of 
inequivalent local minima: 8 corresponding to the 8 sites 
shown in Fig. [I] close to a Gd dopant, and 8 associated 
to the same sites "far" from the dopant, i.e. beyond its 
first oxygen coordination shell. Tab. [T] gives the relative 
energy associated to each site (taken from Ref. [TTj) : in 
the first coordination shell of Gd, only 8 sites among the 
12 can be considered as non-equivalent. Beyond also, 
only the same 8 kinds of sites can be considered as non- 
equivalent with a very good accuracy. In other words, the 
symmetry-breaking caused by the presence of dopants 
can be considered as having no significant influence on 
the energy landscape of the protonic defects. In order to 
distinguish the sites of these two families, we introduce 
another letter, "n" (for a site near the dopant), or "f" 
(for a site far from the dopant). 

To summarize, the 16 kinds of stable positions are la- 
beled by 

• a number (1 or 2) corresponding to the oxygen type 
(apical and equatorial, respectively), 

• a letter ( "a" , "b" , "c" or "d" ) corresponding to the 
0-H direction (shown in figure [l} , 



FIG. 2: Angles between the pseudo-cubic direction and the 
actual O-H bond for the two subcategories of protonic sites. 

• and another letter, "n" for a site near the dopant, 
or "f" for a site far from the dopant. 

In the presence of a dopant, the OH bond might 
slightly deviate from the pseudo-cubic direction: usually 
the proton is expected to bend towards the dopant due to 
the opposite formal charge of the corresponding defects 
(+1 for the protonic defect OHq versus -1 for the dopant 
defect Gd Ce ). But it also depends on the dopant siz e 1 8 1 31 1 . 

In the present case, the proton has indeed a tendency 
to bend slightly towards the dopant, but with a devia- 
tion from the pseudo-cubic direction lower than 10° (see 
Tab. It is possible to divide the eight stable sites 
into two categories: either the proton is able to hop from 
one octahedron to another (a/b-type) or not (c/d-type). 
The c/d-type site shows a noticeable bending (around 5°) 
while the a/b-type are almost perfectly aligned along the 
pseudo-cubic direction. This absence of bending may be 
due to the stabilization of a/b-type sites by an hydro- 
gen bond with the facing oxygen, which is in those cases 
rather close. This hydrogen bond would be dominant 
over the proton-dopant interaction, especially since the 
dopant is much further than for a c/d-type site. 

Note that in perovskites with smaller lattice constant, 
the bending is usually stronger, but also highly dopant- 
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Position 6 near Gd 8 far from Gd 



la 


-0.1° 


0.2° 


lb 


-0.1° 


0.2° 


lc 


5.3° 


0.5° 


Id 


3.5° 


0.5° 


2a 


-0.5° 


0.6° 


2b 


1.6° 


1.2° 


2c 


5.0° 


2.1° 


2d 


8.2° 


4.9° 



TABLE II: Values of the angle described in figure [2] for a 
proton near a dopant, and far from a dopant. 

dependent. Bjorketun et a/P have studied this depen- 
dence in BaZrOa and got a bending angle from 6.9° for 
Gadolinium up to 20.4° for Gallium. An even higher 
bending of around 30° for Scandium, Yttrium or Ytter- 
bium have been found in SrZrO J^U. 



IV. ENERGY BARRIERS 

We have seen that the energy landscape of stable pro- 
tonic sites in Gd-doped BaCe03 is really complex, due 
to the distortions of the Prima structure and the pres- 
ence of dopants. As a result, there are many different 
values for the energy barriers, associated to several diffu- 
sion mechanisms, even by considering the simplified en- 
ergy landscape (with 16 minima) presented previously. 



A. The three different mechanisms: reorientation, 
intra-octahedral and inter-octahedral hopping 

In an ideal cubic perovskite, there are two kinds of pro- 
cesses for theproton motion: reorientation and transfer 
(or hopping}^, to which only two different energy bar- 
riers can be associated, provided the proton is assumed 
to be far from any dopant. In BaZrOa, the reorientation 
(resp. transfer) barrier is 0.14 eV (resp. 0.25 eV), while 
in cubic BaTiOj 311 , it is 0.19 eV (resp. 0.25 eV). In such 
simple systems, each proton in a stable site has four dif- 
ferent possibilities to move: two reorientations and two 
intra-octahedral hopping, the inter-octahedral hopping 
being considered as unlikely (because the oxygen facing 
the OH group is too far). 

However, the existence of tilts of oxygen octahedra, 
very common in perovskite oxides^ having low tolerance 
factor t = rjt+ro makes the inter-octahedral hop- 

V2(rB+ro) 

ping more likely in these strongly distorted structures 
(Fig. [3]), because some inter-octahedral oxygen-oxygen 
distances might be considerably lowered by the anti- 
fcrrodistortive motions of the oxygen atoms. The pro- 
ton may thus jump directly from one octahedron to an- 
other (one inter-octahedral hopping instead of two intra- 
octahedral hoppings), which might result in an increase 
of the macroscopic diffusion coefficients. Tab. |III| em- 
phasizes the link between the tolerance factor t, the per- 




FIG. 3: Possible motions of the proton in the perovskite 
Prima structure, (a) reorientation, (b) intra-octahedral hop- 
ping, (c) inter-octahedral hopping. 



ovskite structure, and the possibility of inter-octahedral 
transfer according to the works mentionned. 

As explained in Sec. |IIIB[ in perovskites with small 
lattice constant (< 4.0 A), the proton in its stable site 
tends to bend towards one oxygen atom of one neighbor- 
ing octahedron instead of being equidistant from both 
neighboring oxygens. In such systems, there are therefore 
twice more stable sites than in perovskites with larger 
lattice constant, so that an additional rotational mech- 
anism might exist, corresponding to the slight reorien- 
tation of OH, bending from the edge of one neighboring 
octahedron to the other. This mechanism was previously 
called "flip"—-' or "bending"^ r "inter-octahedron hop- 
ping'1221 ( Du t "inter-octahedron reorientation" should be 
less confusing, since the bond between O and H is not 
broken during this process). However, the energy bar- 
rier of the flip is usually rather low (< 0.1 eV^, and 
thus most of the time neglected. It can also be seen as 
part of the intra-octahedral transfer mechanism : before 
jumping from one oxygen to another, there is a little re- 
orientation of the proton in order to get an alignment 
O-H...O. The intra-octahedral transfer would thus be a 
two-step mechanism with bending then stretching. 

Tab. |III| illustrates the possible correlation between the 
lattice parameter and the possibility to flip for several 
proton conductor perovskites. Note that some studies 
found a possible inter-octa hedral transfer in small cubic 
perovskite such as SrTiOj^ 14 * 37 ' or even in cubic per- 
ovskites with large lattice cons tant such as BaZrOj^, in 
contradiction with other workiP^'Ml. 
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a (A) 


t 


Structure 


Flip Inter 


4.29 


0.89 


Pnma 


no yes 


4.04 


0.92 


Prima 


no yes 


4.41 


0.94 


Pnma 


no yes 


4.14 


0.95 


Pnma 


no yes 


o.oo 


0.97 


Pnma 


yes yes 


4.25 


1.01 


Pm3m 


no no 


3.91 


1.01 Pm3m (74/ mem) 


yes no 


4.16 


1.03 


Pm3m 


no no 


4.06 


1.07 


Pm3m (Rim) 


no no 



Perovskite 
SrCeO^ 

cazrojsran] 

BaCeOjIl 

Sr z r oj nmin 

CaTiOaPra 

BaZrOsF™ 
SrTiOa™ 
BaSnOjPSl 
BaTiO^ 



TABLE III: Pseudo-cubic lattice constant ao from DFT cal- 
culations (GGA), tolerance factor t calculated from Shannon 
ionic radii, crystal space group of different perovskite oxides, 
and whether flip or inter-octahedral hopping can occur or 
not. For BaTiOs, the high-temperature cubic structure is 
considered, which is the one simulated in Ref. 1351 The cu- 
bic structure is also considered for SrTiOs, rather than the 
low-temperature tetragonal structure. In those two cases, the 
ground state space group is given between parenthesis. 



B. Energy barriers and Minimum Energy Paths 

Using the string method, the Minimum Energy Paths 
joining the various stable sites have been computed, giv- 
ing access to the transition states and thus the energy 
barrier for the corresponding proton motion. These en- 
ergy barriers are provided in Tab. 
riers far from dopants (i.e. from 



Reorientation 



Intra 



Inter 



IV 



Note that the bar- 
to "f ) have been 
also computed in a 80-atom supercell without dopant and 
a +1 charge state (to simulate the protonic defect), com- 
pensated by a uniform charged background. The energy 
barrier values obtained are identical to the ones obtained 
in the doped supercell within 0.01 eV and are presented 
in the Appendix. 

Starting from a given initial position, the possible mo- 
tions for the proton are: two reorientations, two intra- 
octahedral hopping, and possibly one inter-octahcdral 
hopping if the configuration is favorable (which is the 
case for a and b-type positions where the oxygen atom 
facing the proton is close enough). Looking at Tab. IV 
we can notice that barriers between two "near" sites or 
two "far" sites, corresponding to reorientation barriers, 
are very similar (difference within 0.05 eV). This is ex- 
pected as the energy surface of protons bonded to an oxy- 
gen 1st neighbor of a dopant is almost simply shifted by 
0.1 eV compared to that of protons far from the dopant, 
leading to similar energy landscape. However, the case 
of hopping is more complicated since the Coulomb inter- 
action between H and Gd prevents hydrogen from easily 
escaping from the dopant neighborhood. Thus, hopping 
barriers between a "near" site and a "far" site have usu- 
ally a higher value that the ones corresponding to the 
backward motion. 

Fig. [4] illustrates the energy profile for each of the three 
possible kinds of mechanisms (note this is not an exhaus- 
tive list of all possible profiles): Fig. |4^i shows the energy 
profile, as well as the evolution of the O-H distance and 
the angle from the initial O-H direction in the case of a 
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ldn 0.39 2bn 0.13 


2cf 


2df 


0.02 


2bf 


0.17 


2af 


0.36 


2an 0.28 






2df 


2af 


0.20 


2cf 


0.08 


laf 


0.37 


lan 0.34 
















lbf 


0.31 


lbn 0.28 







TABLE IV: Energy barriers (eV) for proton reorientation, 
intra-octahedral hopping ("intra") and inter-octahedral hop- 
ping ("inter"). 



complete turn around an oxygen 01 near the dopant. Us- 
ing the notations of Tab. |IV[ it corresponds to the 4 reori- 
entation mechanisms: lan =>■ lbn len => ldn lan. 
These 4 reorientation barriers have not the same profile 
at all: not only the barrier height can differ by a fac- 
tor 5, but also the angle between two stable sites varies 
from 60° to 120° instead of being set to 90° (case of an 
ideal cubic perovskite). Figs. [IJd and[4j; give similar in- 
formation but for intra-octahedral and inter-octahedral 
hoppings respectively. Both mechanisms seem to occur in 
two steps: first a reorientation, slight for inter-octahedral 
hopping (w 5°) and larger for intra-octahedral hopping 
(w 45°) in order to get O-H-0 aligned, then the jump 
between both oxygen atoms. This reorientation can be 
related to what we mentioned as "flip" in the previous 
section. 



V. DISCUSSION 



Comparison between Gd-doped BaCeC>3 
In-doped CaZrC>3 



and 



The present results on Gd-doped BaCeOs can be 
compar ed with previous values computed in In-doped 
CaZr0 lioi4Q]44| ag hoth materia i s exhibit the same kind 

of structural distortion: BaCeOa and CaZr03 have the 
same perovskite structure with very close Goldschmidt's 
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1b 1c 1d 

Proton position 
(a) Reorientation 



1a 2d 
Proton position 
(b) Intra-octahedral hopping 



1a 1b 
Proton position 
(c) Inter-octahedral hopping 



FIG. 4: Energy profiles and evolution of some geometric quantities along typical Minimum Energy Paths. The angle <j> is 
between the initial and current O-H direction. 





BaCeOs 


CaZrO 3 [l0] 


cubic 


a c (A) 


4.44 


4.06 




a/a c 


1.41 


1.39 


1.41 


b/a c 


1.42 


1.44 


1.41 


c/a c 


2.00 


2.00 


2.00 



A-Q/a e (±<r) 0.71 (±0.21) 0.72 (±0.22) 0.71 
FM5/a c (±<t) 0.51 (±0.00) 0.52 (±0.00) 0.50 
A-O-A (deg) 153.85 144.74 180.00 

B-O-B (deg) 156.45 145.49 180.00 

TABLE V: Structural parameters (lattice parameters, cation- 
oxygen distances and angles) for BaCeOs, CaZr03 and a fic- 
titious cubic perovskite. 



tolerance factor (0.94 and 0.92 respectively) and thus 
have the same orthorhombic structure with Vnma space 
group. However, according to its bigger tolerance factor, 
BaCeC>3 should be slightly less distorted from the cu- 
bic structure and thus inter-octahedral transfer may be 
harder than in CaZrC>3. Tab. |V| confirms that BaCeC>3 
is a bit closer to an ideal cubic structure than CaZrC>3. 



The same tendency is indeed observed with a very large 
range of possible values for energy barriers from a few 
0.01 eV up to nearly 1 eV. For instance, in BCGO, reori- 
entation barriers can take a wide range of different values, 
starting at less than 0.1 eV for barrier between c-type 
and d-type sites up to 0.5 eV for barrier between a-type 
and b-type sites. The same results have been found for 
In-doped CaZrOj^ except for the fact that the largest 
barrier can go up to 0.9 eV. 

The very small barrier between c and d sites might ex- 
plain why position lc is not considered at all in the work 
of Bilic and Gale^ (only 7 different positions instead of 
our 8 positions near a specific B-atom) and 2c near some 
specific oxygen atoms 02. According to Tab.||J lc and 2c 
are much higher in energy than nearby positions, that is 
why the reorientation barriers from c-type site are really 
small. 

In both materials, possible inter-octahedral hoppings 
have a smaller energy barrier than intra-octahedral hop- 
ping. This follows from the ability of any oxygen octahe- 
dron to bend towards another in the orthorhombic Pnma 



structure, so that two facing oxygens (belonging to dif- 
ferent octahedra) can be made very close to each other. 
But each octahedron remains rigid, so that its own oxy- 
gen atoms cannot be made closer to each other (though a 
little distortion during the transfer is observed, in agree- 
ment with previous calculations^). Of course, the inter- 
octahedral hoppings are possible only when the oxygen 
atoms involved are close to each other (this corresponds 
to a/b type within our notation). The c/d type oxygens, 
which are made further from each other as a result of the 
tilting process, are excluded from the inter-octahedral 
motions. 

According to those common tendencies, we can suggest 
that all orthorhombic perovskites behave alike and make 
some assumptions: 

i/ rather low barriers (< 0.2 eV) for inter-octahedral 
hopping depending on the level of distortion (barrier 
is smaller as distortion increases) 

ii/ higher barriers (~ 0.3-0.6 eV) for intra-octahedral 
hopping 

hi/ a wide range of values for reorientation, from less 
than 0.1 eV up to 0.8 eV, depending on the type of 
protonic site. 

Finally, there is a quantitative difference between both 
materials concerning the attractive power of the dopant: 
it seems much harder to escape from Indium in CaZr03 
than from Gd in BaCeC>3. The barrier to escape from In- 
dium is on average three times higher than the backward 
barrier, while in BaCeOs the escaping barrier is higher 
only by 50%. This may be due to the nature of the dopant 
as suggested by Bjorketun et alP, which have shown that 
energy barriers for proton migration near a dopant can 
be strongly dependent of its nature. Therefore Gadolin- 
ium seems to be a good candidate as a dopant since its 
power of attraction is low enough to let the proton escape 
relatively easily. 

B. Rate-limiting events 

The rate-limiting process in such distorted system is 
not so obvious. Contrary to what can be expected, the 
reorientation is not necessarily much faster than the hop- 
ping. Munch and co-workers have found that the proton 
transfer step is indeed rate limiting in BaCeOs but of the 
same order of magnitude as reorientation for SrCeOj^l. 
More precisely, they computed an activation energy for 
rotational diffusion in BaCeOs of 0.07 eV for Oi and 
0.11 eV for O2, close to the values we get for the lowest 
reorientation barriers. In earlier worH^, they found for 
Ba{Ce,Zr,Ti}03 that reorientation happens much faster 
with a time scale of ~ 10~ 12 s, while proton transfer oc- 
curs at a time scale of 10 -9 s. However the three materi- 
als have been studied in their cubic structure, thus pre- 
venting the low-barrier inter-octahedral transfer. Gomez 
and co-workers^2l precise that the rate-limiting process 



in orthorhombic structure is an intra-octahedral trans- 
fer. The fact that most of these studies only focus on the 
cubic structure might explain why the transfer step has 
been thought to be rate-limiting. 



VI. CONCLUSION 

In this work, we have performed density-functional cal- 
culations on fully hydrated Gd-doped barium cerate and 
computed in an exhaustive way the Minimum Energy 
Paths between stable protonic sites close and far from 
the Gd dopant. 

Proton transport in perovskites is usually described as 
a two-step Grotthuss-type diffusion mechanism: a quick 
reorientation, followed by a transfer to another oxygen 5 . 
However, even if this is correct in principle, we have 
found that in Gd-doped BaCe03, the reorientation is 
not necessarily a fast process compared to transfer. In 
this distorted perovskite with orthorhombic Puma space 
group, inter-octahedral hoppings with rather low barri- 
ers ~ 0.2 eV do exist. Also, reorientation mechanisms 
can be very different from one site to another and thus 
take a wide range of possible values from 0.02 eV up to 
0.54 eV. To a lesser extent, the same argument can be 
applied to intra-octahedral hopping for which the energy 
barrier varies between 0.22 and 0.58 eV. 

All these results are qualitatively comparable with a 
previous work focused the orthorhombic perovskite In- 
doped CaZrOs 10 . The low barriers found for inter- 
octahedral hopping in these orthorhombic structures sug- 
gest that protonic diffusion could be much faster in such 
structure than in the cubic one, since an inter-octahedral 
hopping is equivalent to two intra-octahedral transfers 
but with a higher rate. All the barrier values will be ex- 
ploited in Kinetic Monte-Carlo simulations to check the 
actual rate of reorientation versus hopping, and simulate 
proton trajectories on larger space and time scales. 

Finally, gadolinium in barium cerate seems to be in- 
teresting as a dopant as it acts like a shallow trap for 
protons, with rather low escaping barrier (compared to 
indium in calcium zirconate), enabling the proton to dif- 
fuse quite easily. However, other trivalent dopants could 
be tested to check whether they have better properties 
for protonic diffusion. 
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Barrier 


pure 


BaCeOa 


"far" 


BaCeG. 


Reorientation 


— ► 


<- 


— ¥ 




la-lb 


0.54 


0.54 


0.54 


0.54 


lb-lc 


0.33 


0.18 


0.33 


0.18 


lc-ld 


0.06 


0.15 


0.06 


0.15 


ld-la 


0.09 


0.14 


0.08 


0.14 


2a-2b 


0.36 


0.49 


0.36 


0.49 


2b-2c 


0.33 


0.17 


0.33 


0.17 


2c- 2d 


0.03 


0.08 


0.02 


0.08 


2d-2a 


0.20 


0.17 


0.20 


0.17 


Hopping 


— ► 


«— 


— > 


•<— 


la-2d 


0.50 


0.37 


0.50 


0.37 


lb-2d 


0.45 


0.31 


0.45 


0.31 


lc-2b 


0.36 


0.47 


0.36 


0.47 


ld-2b 


0.42 


0.44 


0.42 


0.44 


2a-2c 


0.39 


0.36 


0.39 


0.36 


la-lb 


0.19 


0.20 


0.19 


0.20 


2a-2a 


0.21 




0.21 




2b-2b 


0.16 




0.16 





configuration, have been recomputed using an undoped 
supercell in which the charge of the proton is compen- 
sated by a uniform charged background (jcllium), as fre- 
quently done for the simulation of charged defects. In 
such cases, there are 16 different motions: 8 reorienta- 
tions, 5 intra-octahedral hoppings and 3 inter-octahedral 
hoppings. This corresponds to 30 barrier values. The en- 
ergy barriers obtained using this method are compared to 



the ones obtained using the doped supercell in Tab. [VI 
the values obtained using the two methods are the same 
within 0.01 eV, confirming that a 80-atom supercell is 
large enough to contain a region "close" to the dopant 
and a region "far" from it. The proton "far" from the 
dopant does not feel the influence of Gd atoms, and can 
be considered as in pure BaCe03. Besides, the fact that 
we get the same values in both cases suggests that the 
jellium only induces a systematic shift in total energies, 
but does not affect energy differences. 



TABLE VI: Comparison of barrier values "far" 
dopant in BCGO and in pure charged BaCe03. 



from the 



VIII. APPENDIX: COMPUTATION OF 
ENERGY BARRIERS FAR FROM DOPANTS 
USING A CHARGED SUPERCELL 

The barriers corresponding to motions far from the 
dopant, i.e. from a "f" configuration to another "f" 
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